home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / PowerPC / Dev / PPCRelease / Libm / vax / sqrt.s < prev    next >
Encoding:
Text File  |  1993-09-23  |  4.0 KB  |  124 lines
  1. # Copyright (c) 1985 Regents of the University of California.
  2. # All rights reserved.
  3. #
  4. # Redistribution and use in source and binary forms, with or without
  5. # modification, are permitted provided that the following conditions
  6. # are met:
  7. # 1. Redistributions of source code must retain the above copyright
  8. #    notice, this list of conditions and the following disclaimer.
  9. # 2. Redistributions in binary form must reproduce the above copyright
  10. #    notice, this list of conditions and the following disclaimer in the
  11. #    documentation and/or other materials provided with the distribution.
  12. # 3. All advertising materials mentioning features or use of this software
  13. #    must display the following acknowledgement:
  14. #    This product includes software developed by the University of
  15. #    California, Berkeley and its contributors.
  16. # 4. Neither the name of the University nor the names of its contributors
  17. #    may be used to endorse or promote products derived from this software
  18. #    without specific prior written permission.
  19. #
  20. # THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. # ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. # SUCH DAMAGE.
  31. #
  32. #    @(#)sqrt.s    5.4 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)sqrt.s    1.1 (Berkeley) 8/21/85; 5.4 (ucb.elefunt) 10/9/90"
  38.  
  39. /*
  40.  * double sqrt(arg)   revised August 15,1982
  41.  * double arg;
  42.  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
  43.  * if arg is a reserved operand it is returned as it is
  44.  * W. Kahan's magic square root
  45.  * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82
  46.  *
  47.  * entry points:_d_sqrt        address of double arg is on the stack
  48.  *        _sqrt        double arg is on the stack
  49.  */
  50.     .text
  51.     .align    1
  52.     .globl    _sqrt
  53.     .globl    _d_sqrt
  54.     .globl    libm$dsqrt_r5
  55.     .set    EDOM,33
  56.  
  57. _d_sqrt:
  58.     .word    0x003c          # save r5,r4,r3,r2
  59.     movq    *4(ap),r0
  60.     jmp      dsqrt2
  61. _sqrt:
  62.     .word    0x003c          # save r5,r4,r3,r2
  63.     movq    4(ap),r0
  64. dsqrt2:    bicw3    $0x807f,r0,r2    # check exponent of input
  65.     jeql    noexp        # biased exponent is zero -> 0.0 or reserved
  66.     bsbb    libm$dsqrt_r5
  67. noexp:    ret
  68.  
  69. /* **************************** internal procedure */
  70.  
  71. libm$dsqrt_r5:            # ENTRY POINT FOR cdabs and cdsqrt
  72.                 # returns double square root scaled by
  73.                 # 2^r6
  74.  
  75.     movd    r0,r4
  76.     jleq    nonpos        # argument is not positive
  77.     movzwl    r4,r2
  78.     ashl    $-1,r2,r0
  79.     addw2    $0x203c,r0    # r0 has magic initial approximation
  80. /*
  81.  * Do two steps of Heron's rule
  82.  * ((arg/guess) + guess) / 2 = better guess
  83.  */
  84.     divf3    r0,r4,r2
  85.     addf2    r2,r0
  86.     subw2    $0x80,r0    # divide by two
  87.  
  88.     divf3    r0,r4,r2
  89.     addf2    r2,r0
  90.     subw2    $0x80,r0    # divide by two
  91.  
  92. /* Scale argument and approximation to prevent over/underflow */
  93.  
  94.     bicw3    $0x807f,r4,r1
  95.     subw2    $0x4080,r1        # r1 contains scaling factor
  96.     subw2    r1,r4
  97.     movl    r0,r2
  98.     subw2    r1,r2
  99.  
  100. /* Cubic step
  101.  *
  102.  * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation,
  103.  * a is approximation, and n is the original argument.
  104.  * (let s be scale factor in the following comments)
  105.  */
  106.     clrl    r1
  107.     clrl    r3
  108.     muld2    r0,r2            # r2:r3 = a*a/s
  109.     subd2    r2,r4            # r4:r5 = n/s - a*a/s
  110.     addw2    $0x100,r2        # r2:r3 = 4*a*a/s
  111.     addd2    r4,r2            # r2:r3 = n/s + 3*a*a/s
  112.     muld2    r0,r4            # r4:r5 = a*n/s - a*a*a/s
  113.     divd2    r2,r4            # r4:r5 = a*(n-a*a)/(n+3*a*a)
  114.     addw2    $0x80,r4        # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
  115.     addd2    r4,r0            # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a)
  116.     rsb                # DONE!
  117. nonpos:
  118.     jneq    negarg
  119.     ret            # argument and root are zero
  120. negarg:
  121.     pushl    $EDOM
  122.     calls    $1,_infnan    # generate the reserved op fault
  123.     ret
  124.